function check_convergence(CONSTS, plot_data)

    k0    = CONSTS.k0;
    a     = CONSTS.a;

    m = 5;
    period = pi/(k0*a);
    dots_per_period = 5000;
    num_periods = 40;
    q_vec = (0.0 : period/dots_per_period : num_periods*period)';

    y1_vec = besselmx(double('J'), m+1, k0*a*q_vec, 0).^2;

    y2_vec = (m * (besselmx(double('J'), m, k0*a*q_vec, 0) ./ (k0*a*q_vec))).^2;

    y3_vec = - (m./(k0*a*q_vec)) .* (besselmx(double('J'), m+1, k0*a*q_vec, 0) .* ...
             besselmx(double('J'), m, k0*a*q_vec, 0));

    y4_vec = besselmx(double('J'), m+3, k0*a*q_vec, 0).^2;

    if (plot_data)
        figure; plot( ...
                      q_vec, y1_vec, 'r-', ...
                      q_vec, y2_vec+y3_vec, 'm-', ...
                      q_vec, y3_vec, 'b-', ...
                      q_vec, y4_vec, 'g-'  ...
                      );
        title('func(q)');
        xlabel('q'); ylabel('func'); %ylim([0 8e-13]);
        legend('J_{m+1}^2','(mJ_m/(k_0aq))^2', '-(m/(k_0aq))J_{m+1}J_m', 'J_{m+3}^2', ...
             'location', 'SouthEast');
    end

end